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Testing higher-order Lagrangian perturbation theory 
against numerical simulations - 

1. Pancake models 

by 

T. Buchert, A.L. Melott, A.G. Weifi 


Summary: We present results showing an improvement of the accuracy of perturba- 
tion theory as applied to cosmological structure formation for a useful range of quasi- 
linear scales. The Lagrangian theory of gravitational instability of an Einstein-de Sitter 
dust cosmogony investigated and solved up to the third order in the series of papers 
by Buchert (1989, 1992, 1993a), Buchert & Ehlers (1993), Buchert (1993b), Ehlers &: 
Buchert (1993), is compared with numerical simulations. In this paper we study the 
dynamics of pancake models as a first step. In previous work (Coles et ah 1993, Melott 
et al. 1993, Melott 1993) the accuracy of several analytical approximations for the 
modeling of large-scale structure in the mildly non-linear regime was analyzed in the 
same way, allowing for direct comparison of the accuracy of various approximations. In 
particular, the “Zel’dovich approximation” (ZePdovich 1970, 1973, hereafter ZA) as a 
subclass of the first-order Lagrangian perturbation solutions was found to provide an 
excellent approximation to the density field in the mildly non-linear regime (i.e. up to a 
linear r.m.s. density contrast of a ^ 2). The performance of ZA in hierarchical cluster- 
ing models can be greatly improved by truncating the initial power spectrum (smoothing 
the initial data). We here explore whether this approximation can be further improved 
with higher-order corrections in the displacement mapping from homogeneity. We study 
a single pancake model (truncated power-spectrum with power-index n — —1) using 
cross-correlation statistics employed in previous work. 

We found that for all statistical methods used the higher-order corrections improve 
the results obtained for the first-order solution up to the stage when <r(linear theory) 
% 1. While this improvement can be seen for all spatial scales, later stages retain 
this feature only above a certain scale which is increasing with time. However, third- 
order is not much improvement over second-order at any stage. The total breakdown 
of the perturbation approach is observed at the stage, where cr(linear theory)^ 2, 
which corresponds to the onset of hierarchical clustering. This success is found at 
a considerably higher non-linearity than is usual for perturbation theory. Whether a 
truncation of the initial power-spectrum in hierarchical models retains this improvement 
will be analyzed in a forthcoming work. 
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1. Introduction 


In most recent research in cosmology, the gravitational instability picture of the standard 
Friedman- Lemait re cosmogonies has been used to describe the earliest stages of the forma- 
tion of inhomogeneities in the Universe (see, e.g, Peebles 1980 and ref. therein). The first 
quantitative exploration used linear perturbation theory to make predictions from theory. 
Although there are serious restrictions of applicability of this perturbation approach (e.g., 

the density contrast 6 = p ~ pii , where pu is the homogeneous density of the background 

cosmogony, should be much less than 1), the solutions are still applied to normalize inho- 
mogeneous cosmological models, i.e., to determine the amplitude of the initial fluctuation 
field in comparison with the density contrast we observe today. We now know that first- 
order solutions of this perturbation approach cannot be used as an approximation of the 
present day density field. Only, if we smooth the field with an extremely large smoothing 
length of the order of 100 h -1 Mpc (Coles et al. 1993), a comparison with non-linearly 
evolved density fields can be meaningful. This perturbation approach is Eulerian, i.e., the 
perturbations are evaluated as functions of Eulerian coordinates. 

In the 70’s Zel’dovich realized the limitations of this approach and proposed an ex- 
trapolation of the Eulerian linear solutions into the mildly non-linear regime by employing 
the Lagrangian picture of continuum mechanics (Zel’dovich 1970, 1973). In his model 
the motion of the continuum is similar to the purely kinematical motion of a system of 
particles which move under inertia. The model can be mapped to this system by a suit- 
able transformation of dependent and independent variables (Shandarin & Zel’dovich 1984, 
1989). The success of Zel’dovich’s extrapolation for the description of typical patterns of 
the large-scale structure has been widely recognized, although not generally understood. 
It can be traced back to the fact that the Lagrangian description is based on following 
the trajectories of fluid elements, while the density itself does not appear as a dynamical 
variable in this picture and can be integrated exactly along the trajectories of a given 
solution. Consequently, the density is not limited in a Lagrangian perturbation approach, 
only the gradient of the displacement field has to be small compared with the expansion 
factor of the homogeneous background (Buchert 1992, Ehlers & Buchert 1993). Moreover, 
Zel’dovich’s model can be systematically derived by solving the Lagrangian evolution equa- 
tions to the first-order (Buchert 1992). It also provides a class of exact three-dimensional 
solutions in the case where the collapse is maximally anisotropic (i.e., in the case where lo- 
cally two of three eigenvalues of the fluid’s deformation tensor vanish) (Buchert 1989). The 
plane-symmetric case is contained in a subclass of these solutions; first-order Lagrangian 
perturbation solutions are the general solutions of the full system in this case. 

Lagrangian perturbation solutions are now available for various backgrounds and up 
to the third order in the deviations from homogeneity (see Buchert 1993b and references 
therein). Parallel efforts to investigate Lagrangian perturbation solutions are made by 
Moutarde et al. (1991), Bouchet et al. (1992), Lachieze-Rey (1993), Gramann (1993), 
Giavalisco et al. (1993), Croudace et al. (1993), see also Bertschinger & Jain (1993). 

Coles, Melott and Shandarin (1993) conducted a series of tests of analytical approxi- 
mations by cross- correlating them with N-body simulations. They tested the Eulerian lin- 
ear perturbation solution and the “Zel’dovich approximation” (hereafter ZA) as a subclass 
of the irrotational Lagrangian linear solution. They found that the latter was considerably 
more successful than the former. Considerable further improvement is made if the ZA is 
applied not to the full power spectrum on the whole range of spatial scales, but only to a 
truncated body of the spectrum. The truncation removes unwanted non-linearities, which 
are evolved beyond the range of applicability of the model (Melott and Shandarin 1990, 
Beacom et al. 1991, Kofman et al. 1992, Coles et al. 1993). The shape of high-frequency 
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filters imposed on the power-spectrum was found to be important (Melott et ai 1993). 
Most successful in respect of improvement of the cross-correlation statistics was the use of 
Gaussian filters, but if one wishes to work backwards from evolved to initial state, a sharp 
truncation (step function) in k is preferable (Melott 1983). 

Since closer study of analytical approximations reveals both limitations and powerful 
properties, it is worthwhile to study generalizations of Zel’dovich’s mapping to improve the 
cross-correlation with numerical simulations. In this way the modeling of large-scale struc- 
ture is effective and can be used to provide initial and boundary conditions for the evolution 
of small-scale structure within the large-scale environment, which could be modelled, e.g., 
by fully hydrodynamical simulations. Although the applicability of the analytical models 
under study is not limited to a particular type of spectrum (the “Zel’dovich approxima- 
tion” was originally used for modeling the standard “pancake picture”, see Shandarin &: 
Zel’dovich 1989), we here do concentrate our comparison on “pancake models”, i.e., models 
which do not involve small-scale power in the initial data. We do this as a first step to 
properly understand the limitations of higher-order corrections to ZeFdovich’s mapping. 
Also, pancake models can be understood as generic archetype of hierarchical models which 
involve pancaking on all spatial scales (Kofman et ai 1992). Because the models we study 
here already have a truncated initial power spectrum, we do not apply additional trunca- 
tion as in Coles et ai (1993) or Melott, Pellman and Shandarin (1993). In a second step, 
we shall later analyze various power spectra which are evolved deeply into the non-linear 
regime in order to probe the performance of the Lagrangian approximations in the case of 
hierarchical models (paper II in preparation ). 

Further work on other approximations is also in preparation: the frozen-flow approx- 
imation (Matarrese et ai 1992) and the adhesion approximation (Kofman et ai 1992 and 
references therein) will be given a similar treatment. 

Another advantage of analytical models is the possibility of high-spatial resolution of 
the density field. Analytical mappings of a continuum offer in principle infinite resolution 
in contrast to numerical N-body simulations. We shall explain that models with suffi- 
ciently smooth initial data like pancake models can be resolved easily to, e.g., an effective 
resolution of 2048 3 particles. However, as we shall see later, the accuracy of this effort 
would be limited to early stages of non-linearity. 


2. Numerical realization of pancake models 

We specify initial data in terms of a power spectrum V(k) (as a function of comoving 
wavenumber k = \k\) of Gaussian density perturbations of the form: 

V(k)=< |^| 2 >oc|fc|" i (1) 

where 8^ is the discrete spatial Fourier transform of the density contrast 6 . We shall take 
n = — 1, and truncate the power spectrum at a characteristic frequency k c . In this way we 
will probe the epoch of first pancaking. We also define the non-linearity scale k n i(t): 

fibril ( t) 

a 2 {t) / d 3 kV(k) = 1 , (2) 

Jo 
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where k„i(t) is decreasing with time as successively larger scales enter the non-linear regime; 
a(f) is the scale factor of the homogeneous background: for the initial data we set k c = 4k f, 
and a(f*) := 1 where kj is the fundamental mode of the simulation box. We emphasize that 
we give initial data early (at the r.m.s. density contrast of cri(k c ) = 0.01) to guarantee an 
objective modeling of the collapse of first pancakes in the model. One of us (Melott 1985, 
1987) has emphasized the importance of beginning numerical simulations at a highly linear 
stage to avoid an artificial delay of the collapse time of first objects. This delay is inherent 
in Zel’dovich’s displacement mapping which is used to initialize the numerical simulation 
(compare Blanchard et al. (1993) for a discussion of this problem). From several first 
studies of higher-order perturbation solutions we know that the collapse is significantly 
accelerated by the higher-order corrections even in a generic field (Buchert et al. 1993). 
An objective test of this feature can only be conducted by giving the initial data early. 
The results of making too late a start with too high an amplitude are not detectable in the 
two-point correlation function since this is mostly determined by long waves just going 
non-linear (Beacom et al. 1991; Melott & Shandarin 1993). 

The evolution of inhomogeneities is modelled in a flat Einstein-de Sitter background 
( a(t ) = (jf:) 2 / 3 ). We evolve the field with an enhanced PM (particle-mesh) method (Melott 

1986) using 128 3 particles each on a comoving 128 3 mesh with periodic boundary condi- 
tions. This method makes the code resolution-equivalent to a traditional PM code with 
128 3 particles on a 256 3 mesh (see Weinberg et al. (1993)). 

In what follows we will specify the stages we study by the Eulerian linear theory 
amplitude which is just proportional to the scale factor in these fi = 1 models. (The r.m.s. 
density contrast 0*2 denotes ^(linear theory) in the numerical simulation throughout this 
paper.) Of course, the real density contrast is larger due to mode coupling effects. 

From this simulation we extract six stages, the first at the epoch when pancake forma- 
tion is just beginning (02 = 0.25), the second at the moment when the first pancakes have 
collapsed ( 0 ^ = 0.5), the third at a stage when the cellular connection of different pan- 
cakes is established (corresponding to the contemporary epoch according to the standard 
normalization) (02 = 1)? the fourth and fifth at more advanced stages, where secondary 
shell-crossings have developed (02 — 1.2 and cr 2 — 1.5), and the last at a stage when the 
non-linearity scale has dropped to k n i = 2k t ( 0*2 — 2). This stage corresponds to the onset 
of hierarchical clustering due to merging of pancakes. 


3. The Lagrangian perturbation solutions 


In what follows we make use of a restriction of the initial data commonly used in realiza- 
tions of large-scale structure models as well as in the numerical realization quoted above. 

We require that, initially, the peculiar- velocity u(X,ti) be proportional to the peculiar- 
acceleration w(X ,ti): 

u(Xyti) = w(X,ti)ti , (3) 

where we have defined the fields as usual (compare Peebles 1980, Buchert 1992). This 
restriction has proved to be appropriate for the purpose of modeling large-scale structure 
since, for irrotational flows, the peculiar- velocity field tends to be parallel to the gravi- 
tational peculiar-acceleration after some time. The reason for this tendency is related to 
the existence of growing and decaying perturbations in the linear regime, the growing part 
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supports the tendency to parallelity. The assumption of irrotationality should be adequate 
down to the non-linearity scale. 

Henceforth, a comma denotes derivative with respect to Lagrangian coordinates, and 
A 0 := Vq, where Vo denotes the nabla operator with respect to Lagrangian coordinates. 
We shall use in the following the initial peculiar- velocity potential S defined as u(X,ti) 

Vo S(X). The initial peculiar-gravitational potential <f>, w(X ,ti) =: — Vo<^(A’’) is related 
to it as S = —<f>U (eq. (3)). Further, we introduce the symbols /(£,*,*) = tr(S t i ik ) = A 0 <S, 
II{S, ijk ) = |[(tr(5 iifc )) 2 - M(S, i)fc ) 2 )] and III(S^ k ) = det (£,<,*) which denote the three 
principal scalar invariants of the tensor (5^*). 

We introduce the field of trajectories x = /(A^,/), where x denote Eulerian and X La- 
grangian coordinates. We express the solutions in a comoving reference system q = F(X , t), 
where q = x/a(t) denote comoving Eulerian coordinates. 

With a superposition ansatz for Lagrangian perturbations of an Einstein-de Sitter 

background the following family of trajectories q = F(X,a) as irrotational solution of the 
Euler-Poisson system up to the third order in the perturbations from homogeneity has 

been given by one of us (Buchert 1993b). The general set of initial data (<f>(X) i S(X)) is 
only restricted according to S = — <f>ti (eq. (3)): 

F = X + q,(a)V 0 S^(X) + q 2 (a)V 0 S^ 2 \X) 

+ q 3 (a) Vo 5 (3o) (X) + 93(a) V 0 S (36) (X) - q c 3 {a) V 0 x 5 (3c) (X) , (4) 

with the time-dependent coefficients expressed in terms of the expansion function a(i) = 

(^) 2/3 : 



(4a) 

(46) 

(4c) 

(4d) 

(4e) 


The perturbation potentials have to be constructed by solving iteratively the following set 
of 7 Poisson equations: 


A 0 S^ = I(S ii<k )ti , 

(4/) 

A 0 S< 2 > = 211 ( 5 %) , 

(4<?) 

A 0 «S (3a) = 3 III{sV k ) , 

(4 h) 
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(40 


A 0 S< 36 > 


apw.syUe) 

^ €abc d(X u X 2 ,X 3 ) 

a } b t c 


(A„S»% = <„ (i 


e(x„x„x,) 


(4j, fc, /) 


(i,j,k = 1,2,3 with cyclic ordering). 


We realize the solution by first solving Poisson’s equation for S via FFT (Fast Fourier 
Transform) from the initial density contrast 6 generated as initial data for the numerical 
simulation (Section 2). In a flat background model we have: 

A,S = -~S. (5) 

We first note that the first-order part of the mapping (4) is equivalent to Zel’dovich’s 
approximation, if we put = Sti . The only formal difference is that we start at the 

initial time on an undeformed grid due to the assumption a(t{) = 1: F(X,ti) = X . This 
solution of the first Poisson equation (4f ) as well as the entire solution (4) is unique provided 
we impose periodic boundary conditions on S and fix some gauge conditions accordingly 
(Ehlers & Buchert 1993). With this setting the source expressions in (4f-h) are computable 
in terms of the initial data S . We solve all derivatives and all Poisson equations in k — space 
by using FFT, the four Poisson equations (4i-l) have to be solved iteratively by using the 
second-order perturbation potential calculated first from (4g). 

The first-order part of the solution (4) covers the kinematical aspect of the collapse 
process, whereas the second-order part covers essential aspects of the tidal action of the 
gravitational field (Buchert &; Ehlers 1993). At the third order, there are also interaction 
terms among the perturbation potentials at first and second order, which induce vorticity in 
Lagrangian space. Although the solution (4) is rotational in Lagrangian space, the Eulerian 
fields remain irrotational until the first shell-crossings. The development of caustics will 
result in multi-stream flow and will introduce vorticity also in Eulerian space (see, e.g., 
Chernin 1993). Note that the Eulerian representation of the equations breaks down at 
caustics, while the Lagrangian representation of the fluid’s motion in terms of trajectories 
is still regular and formally allows to follow the motion of the fluid across caustics which 

we shall do. However, we emphasize that the Lagrangian solution / represents solutions of 
the Euler-Poisson system only as long as the mapping to Eulerian space is single- valued; In 

multi-stream regions a formed extrapolation of / across caustics neglects the gravitational 
interaction of the streams since, in these regions, the gravitational field-strength acting on 
each particle is the sum of the field-strengths of the streams (compare Ehlers & Buchert 
1993). 


In order to understand that the higher-order approximation schemes are better than 
the first-order scheme until shell-crossing (regardless of the initial conditions and the scale 
of the fluctuations), the following general equation can be used: Consider the density 
contrast A := (p — Qh)/ Q, — oc < A < 1, which is more adapted to the non-linear situation 
than the conventional definition 6 = (g — Qh)I 6h = A/(l — A), defined in Eulerian 
perturbation theory. For this the following fully non-linear evolution equation derived 
from the continuity equation and Poisson’s equation has been found (Buchert 1992; see 
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also Peebles 1987, Buchert 1989, Bertschinger & Jain 1993): 


(ба) 

(бб) 


A = (A - 1)1 , 

A + 2- A - 4nGg H A = (A - 1)211 , 
a 


where here I and II denote the first and second principal scalar invariants, respectively, 
of the peculiar- velocity tensor gradient (uij) (contrary to the notation before, a comma 
here denotes derivative with respect to comoving Eulerian coordinates q, a dot denotes 
Lagrangian time-derivative ^ := -j^\ q + ^ • V g ): 


a X u ‘-‘ ’ 2 a 2 (<5>M) 2 Eta)) * 


(6c) 


This equation is exact. For the approximation schemes we can read off the following: The 
first-order model solves equation (6) with II — > 0, whereas the second-order model takes 
the second invariant into account lor small displacements from the first-order trajectories. 
Thus, the second-order theory improves upon first-order until equation (6) breaks down 
at shell-crossing. 


We realize the solution (4) at each order by calculating 128 3 trajectories. Note that, 
although we have 7 Poisson equations to solve, the model (4) is very time-efficient: On 
a CONVEX C3 the third-order model takes about 20 minutes to realize the density dis- 
tribution for each stage at this resolution, while the first-order model takes 2 and the 
second-order model 3 minutes. On a CRAY YMP/464 the corresponding CPU times are 
smaller by a factor of 5. The N-body simulation run here takes 138 integration steps each 
of which uses about 2 minutes on a CRAY 2, i.e. in total 5.4 hours to get to the final 
stage. The N-body code is limited by the speed of the (rather poor) FFT on a CRAY 2. 
Since this was written, we have moved our N-body computing to a CONVEX C3. 

A higher resolution can be easily obtained for the analytical model, provided the 
smallest wavelength in the model (corresponding to k c = 4k f here) is calculated to suf- 
ficient accuracy. For example, if we resolve the box with length L with 128 calculated 
density values along one spatial direction, the smallest wavelength Lf 4 will be calculated 
at 32 points; 32 points will be sufficient to map the coherence scale at high accuracy. 
Interpolating the calculated grid points, e.g., 16 times, we obtain an effective resolution 
in the box of 2048 3 particles. This interpolation method gives highly accurate results for 
the final distribution, since the analytical mapping contains all non-linear information; it 
maps the interpolated grid points as if it were calculated ones. For a detailed discussion 
and illustration of this method see (Buchert & Bartelmann 1991). 

Figures la to Id show a comparison of slices of the N-body model and the first- through 
third-order approximate schemes by visual appearance. Often very high-order information 
which is difficult to quantify can be rapidly gained from visual presentations. 
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4. Cross-correlation statistics 


4.1. Cross-correlation coefficient 


As in Coles et al. (1993) and Melott et al. (1993) we use here the usual cross- correlation 
coefficient S to compare the resulting density fields: 


S := 


< (M2) > 

<TlCT 2 


(7) 


where 8t,l = 1,2 represent the density contrasts in the analytical and the numerical 
approximations, respectively, at = \/< 6j > — < 8t > 2 is the standard deviation in a 
Gaussian random field; averages < ... > are taken over the entire distribution. 

The density in the analytical models is calculated by collecting trajectories of the 
Lagrangian perturbation solutions at the different orders into a 128 3 pixel grid with the 
same mesh method (CIC binning) as in the N-body simulation. 

The correlation coefficient obeys the bound |S| < 1; 5 = 1 implies that 8\ — C 8 2 , with 
C constant; C close to 1 implies, of course, better agreement between the approximation 
and the numerical simulation. 

Following Coles et al. (1993) we will compare the approximation with the simulation 
with varying amounts of Gaussian smoothing of each so that we can follow accuracy at 
various lengthscales, and allow for the possibility of results which are essentially correct 
but for small displacements. We will plot S as a function of a 2 , the r.m.s. density contrast 
in the numerical simulation. The cross-correlation coefficient S is calculated down to the 
resolution scale of a 128 3 pixellization of the density field. 

There will be basic differences from the procedure of Coles et al. . No truncation of 
the spectrum will be considered here because in pancake models the initial conditions are 
already truncated. We will not be able to compare these results with those since we are 
dealing with a different class of models. Instead, we will compare the accuracy of first— 
, second-, and third-order Lagrangian approximations with the simulation, and thereby 
gain information on the benefits of going to higher-order Lagrangian solutions within 
the context of pancake models. In paper II on hierarchical models we will determine to 
what extent these results carry over to hierarchical models, and we will find out whether 
substantial improvements over the optimally truncated first-order approach (Melott et al. 
1993) is made. 

In order to couple our cross-correlations to the physical scale, we show the actual r.m.s. 
density fluctuations for all approximations at as a function of the Gaussian smoothing 
scale in grid units in Figure 2 for each of our stages. The Gaussian smoothing is done 

by convolution with e~ R / 2R ° . The stages correspond to (Eulerian) Unear theory a 2 — 
0.25,0.5,1.0,1.2,1.5, and 2.0 as discussed before. Rq is given in grid units; remember 
that the box is 128 grid units in size, and the smaUest wavelength perturbation in the 
initial data was 32 grid units. Since it is known that these Lagrangian schemes give rather 
accurate locations for pancakes, we are more interested in studying details of structure. 
Therefore most of our information will be relevant for Rg < 32 grid units. 

In Figures 3a to 3f we show our most important basic result, the cross-correlation 
coefficient of the various orders of approximation with the N-body simulation at various 
stages. 
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4.2. Power spectrum analysis 

Since Lagrangian perturbation solutions include substantial non-linearities (see, e.g., 
Buchert 1992), we expect that they do not preserve the initial power spectrum of the 
fluctuation field as the Eulerian linear solution would do. Therefore, it has meaning to 
test the agreement with the N-body simulation. For the spectrum considered in the present 
paper (n = —1) Melott et al. (1993) found that the power spectrum of the truncated ZA 
with that of the N-body simulation is substantially underestimated. In this respect the 
spectrum with index n = — 1 provides an extreme case. We, here, analyze whether the 
higher-order corrections can improve on that, and present the results of the power spectrum 
analysis in Figures 4a to 4f. 

4.3. Relative phase-errors 

The non-linear evolution of the fluctuation field does not preserve the initial power spec- 
trum as well as the phase information, which is not contained in the power spectrum. There 
will be a substantial change in the phases due to non-linear phase-locking effects, which 
are already described by the first-order Lagrangian perturbation solutions (see Melott et 
al. 1993). Also here we expect an improvement due to higher-order corrections, since they 
contain essential effects of tidal interactions which should contribute to the change of the 
phases. 

We calculate the relative phase error as follows: The Fourier coefficients of the initial 
fluctuation field 8 g = contain information about the amplitude |£g| and the phase 

angle We measure the angle 0 = <pi — <f> 2 , where <j > 2 are the phases in the N-body 
simulation, and <f>i in the analytical approximations. We present the results on the relative 
phase errors in terms of cos(0) as a function of k (calculated in spherical shells in k — space). 
The opposite assignment 6 = <f > 2 — <f>\ yields equivalent results since cos(0) is symmetric 
around 0 = 0. Perfect agreement between the N-body result and the analytical scheme 
implies cos(0) = 1, anti-correlated phases have cos(0) = —1, and for randomized phases 
cos(0) would average to 0. We present the result in Figures 5a to 5f. 


4.4. Density distribution functions 

The cross-correlation coefficient (4.1) measures differences in position up to a constant 
ratio C of the compared densities (compare Section 2). In order to probe the differences 
of the various approximations with respect to the actually predicted densities on different 
spatial scales, we calculate the density distribution function for each approximation and 
depict them in Figures 6a to 6d in terms of the number of cells N found with mass density 
p (in units of the mean) as a function of scale. 


5. Discussion of the results 


In Figures la to Id we show slices of the density fields corresponding to the stages 2,3,4 
and 5, while the statistical analysis is done for all six stages. 

At the first stage (<T 2 = 0.25) pancake formation is just beginning. Here, the higher- 
order schemes predict structures which are almost identical to the N-body result, whereas 
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the first-order scheme appears to delay the collapse time as was expected. The cross- 
correlation coefficient as well as the power spectrum measures this feature clearly (Figs. 
3a, 4a); the phase information is represented better at higher orders especially on small 
scales (Fig. 5a). Also the density distribution functions confirm this behavior (Fig. 6a). 
The number of high-density cells is underestimated by all schemes. 

Stage 2 ( <r 2 = 0.5, Figs, la) emphasizes this: Here, the higher-order schemes significantly 
improve on the first-order “Zel’dovich approximation”; the cross-correlation coefficient 
(Fig. 3b) is as high as S = 0.90 for the second-order and only slightly higher for the third- 
order scheme, while the first-order coefficient has 5 = 0.85 on the scale where a 2 = 1. A 
similar conclusion can be drawn from the power spectrum analysis (Fig. 4b). However, 
the power on small scales is already underestimated by all schemes, but represented better 
at higher orders. Interestingly, the phase information contained in the first-order approx- 
imation shows randomized relative phase-errors, while the higher-order schemes show a 
strong tendency to anti-correlated phases at this stage (Fig. 5b). This point would be- 
come clearer by analyzing the velocity fields. The better agreement with respect to the 
cross-correlation and the power spectrum is also a consequence of a better prediction of 
the collapse time in the higher-order solutions. The density distribution functions (Fig. 
6b) show good agreement with the N-body distribution for the second- and third-order 
models. The first-order model still suffers from its properties mentioned before at stage 1. 

At stage 3 (corresponding to the present stage according to the standard normalization, 
<r 2 = 1) we can already see slight differences in the density fields (Figs, lb): The biggest dif- 
ference is apparent in the lower-right corner where the most evolved structure contained in 
this slice is seen. We appreciate that in the higher-order solutions filaments stay more com- 
pact like in the N-body slice in contrast to the first-order solution. The cross- correlation 
coefficient (Fig. 3c) shows an almost constant improvement of the second-order upon the 
first-order approximation down to the smallest scales, while the third-order scheme be- 
comes worse than the second-order result on scales a 2 > 2, but still is better than that 
of first-order up to <r 2 % 4. The power spectrum (Fig. 4c) remains better represented in 
the higher-order solutions. Starting from this stage, the relative phase-errors show a ten- 
dency to randomization for all orders of the perturbation solutions, while the scale where 
randomization occurs increases with time (Figs. 5c-f). The density distribution functions 
(Fig. 6c) show excellent agreement with the numerical result but begin to overestimate the 
number of cells with moderately high density, a feature which is due to the overestimated 
size of the pancakes. 

Stage 4 (<t 2 = 1.2) is the most interesting stage, where all orders of approximation schemes 
provide equally good results. The agreement of the density fields with the N-body simula- 
tion is still reasonable in all schemes. The most evolved structure shows compact internal 
density peaks inside the first-order pancakes (Figs. 1c). The cross- correlation coefficients 
(Fig. 3d) are still high, but only slightly higher in the second- and third-order approxima- 
tions on large scales corresponding to a 2 < 3 (second-order), cr 2 < 2.5 (third-order). On 
smaller scales the second-order coefficient stays close to the first-order coefficient, while 
the third-order scheme starts to fall off more rapidly below the cross-correlation of the 
first-order scheme. The power spectra (Fig. 4d) are almost identical for all schemes, a 
property which roughly remains in the following stages. Here, the density distribution 
functions (Fig. 6d) mirror an amplification of the effect mentioned for stage 3. 

The next stage 5 ( <r 2 = 1.5) we can identify with the stage where the break-down of the 
Lagrangian perturbation solutions has started. Figs. Id show differences in the thickness 
of pancakes, which begin to grow more rapidly in the perturbation solutions than in the 
N-body simulation. The density distribution functions (Fig. 6e) overestimate the number 
of low-density cells, while the number of high-density cells is underestimated. The scale 
above which an improvement of the higher-order schemes can be detected shifts to larger 


10 



scales, here the scale is roughly <r 2 « 1.5 (Fig. 3e). Furthermore, the improvement is now 
very small. The power spectrum of higher order schemes is now than worse than first-order 
on small scales (Fig. 4e). 

At stage 6 (a 2 = 2.0) such a scale no longer exists, and we cannot appreciate any improve- 
ment of the higher-order schemes upon first-order. The break-down of the perturbation 
sequence is reached. To summarize, stage 3 ( a 2 — 1.0) is the last stage at which the im- 
provement of second- and third-order schemes over first is large and positive. Increasingly 
in stage 4 and later differences become insignificant. There is substantial benefit at early 
stages from second-order, but there is never any important additional improvement by 
going to the third-order. 


6. Conclusions 


We have analyzed a time sequence of stages in a pancake model simulated numerically 
and compared with analytical Lagrangian perturbation solutions at various orders. We 
found that until the stage when ^(linear theory )= 1 (corresponding to the present epoch 
according to the standard normalization of large-scale structure models), the higher-order 
Lagrangian solutions clearly improve upon the first-order “Zel’dovich approximation” down 
to the resolution scale of the simulation. At later stages any small improvement can 
only be appreciated above a scale which increases with time until the stage ^(linear 
theory )% 1.5. Later stages correspond to the onset of hierarchical clustering, which are 
beyona the reach of the perturbation solutions at any order in this model. Also, third- 
order is computationally different and does not improve things much beyond second— order 
at any stage. 

A typical feature of this break-down is the behavior of the higher-order schemes, 
which both become worse than the first-order scheme at later stages. This is due to the 
fact that in a perturbation approach the higher-order time coefficients grow more rapidly 
than the first-order time coefficient (here: q\ oc a, q 2 ex a 2 , oc a 3 ). Consequently, the 
higher-order corrections blow up at and after the time coefficient functions are of the same 
order in all schemes which is the case roughly at the stage ^(linear theory)— 1.5. The 
limit (or break-down) scale above which an improvement of the higher-order schemes 
upon the “Zel’dovich approximation” can be detected at about a 2 ( k 5 ) ~ 2 corresponds to 
Rg ~ 5 and shifts to larger scales at later stages. After this stage we probe the onset of 
hierarchical clustering, for which the perturbation solutions are not meant for in the first 
place. Whether a truncation of the power spectrum in hierarchical models can retain the 
improvements detected here, will be analyzed in paper II. We also think that the break- 
down of the approximations must not be attributed to the lack of higher-order corrections. 
It might well be that the neglection of the self-gravitating interaction of streams in multi- 
stream regions by simply extrapolating the trajectories across caustics must be considered 
a major source of this break-down phenomenon. 

We emphasize that both the N-body simulation and the Lagrangian perturbation 
solutions are assumed to be approximations to the unknown exact solution. 

Another result of this paper is the relative importance of terms in the higher-order 
schemes, if we are concerned with the large-scale performance of the solutions in pancake 
models: We found that the vector perturbation potential can be safely neglected 

in this case. It will only be relevant if we resolve internal structures of pancakes. The 
same can be said for the third-order correction as a whole which does not improve upon 
the second-order scheme on large scales except with respect to the density distribution 
function at early stages. However, internal structures of pancakes are affected due to the 
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prediction of a third shell-crossing inside pancakes; also the collapse occurs substantially 
earlier than in the second-order model if calculated at high resolution (see Buchert et ai 
1993). The higher core density in clusters as a result of this does not yield a noticeable 
improvement in the cross-correlation with the N-body simulation in the present case except 
a slight improvement of the cross-correlation of the third-order model on small scales at 
later stages. We also found that the realization of the third-order model is extremely 
sensitive to numerical errors by iteratively solving Poisson equations. A direct analytical 
calculation by explicitly solving the Poisson equations for a given model is more reliable 
but extremely CPU time extensive for a generic fluctuation field (see: Buchert et al. 1993). 
The second-order scheme which is much easier and faster to realize provides the main effect 
of improvement upon the “Zel’dovich approximation”. It also provides a good compromise 
between improvement at early stages, while the surprisingly good performance of the first- 
order scheme on small scales at later stages is essentially preserved. Here we emphasize that 
the first impression of moderate improvement of second upon first order is misleading, since 
we compare with an approximation which already shows excellent agreement with the N- 
body result: e.g., the cross-correlation coefficient is already high in the first-order scheme. 
The latest stages analyzed here are in a regime where we are moving into hierarchical 
clustering. A removal of unwanted non-linearities by using, e.g., a Gaussian smoothing 
of the initial data will definitively change this picture (MPS). We will examine this when 
we study hierarchical models in paper II. We finally note that second-order solutions are 
available for a larger class of initial data than considered in this paper (Buchert & Ehlers 


1993b 

To summarize, without reservation we recommend the use of second-order perturba- 
tion theory up to the stage cr(linear theory )= 1, as a definite and useful improvement upon 
first-order (the “Zeldovich approximation”) in pancake models. 


The numerical code to realize Lagrangian perturbation solutions is available via e-mail: 
TOB @ ibma.ipp-garching.mpg.de . 
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Figure Captions 


Figure 1: A thin slice (thickness L/128) of the density field is shown for the numerical (upper 
left), the first-order (upper-right), the second-order (lower left), and the third-order (lower right) 
approximations for the evolution stages 2 (<r 2 = 0.5, Fig. la), 3 (02 = 1.0, Fig. lb), 4 (<r 2 = 1.2, 
Fig. lc), and 5 (<r 2 = 1.5, Fig. Id). The grey-scale is linear, and the maximum density contrasts 
are chosen to be 8 in Fig. la, and 25 in Figs. lb,c,d. (The slight stripes in low-density regions are 
artifacts resulting from the interaction of the distorted grid of the particles with the pixelization.) 

Figure 2: The standard deviations of the density contrast as a function of smoothing scale Rq 
in the first-order (dotted), the second-order (dashed), the third-order (dashed-dotted), and the 
numerical (full line) approximations, for the evolution stages 1,2, 3, 4, 5, 6 (Figs. 2a,b,c,d,e,f). Rq 
is given in grid units; the shortest wave in the initial data was 32 grid units. 

Figure 3: The cross-correlation coefficient 5 as a function of the standard deviation er 2 for the 
evolution stages 1,2, 3, 4, 5, 6 (Figs. 3a,b,c,d,e,f). The cross-correlation of the N-body with first- 
order perturbation theory is shown as a dotted line; with second-order a dashed line; and with 
third-order a dashed-dotted line. 

Figure 4: The power spectra of the N-body simulation (solid line) compared with the first- 
dotted), second- (dashed), and third-order (dashed-dotted) approximation schemes for the stages 
1,2, 3, 4, 5, 6 (Figs. 4a,b,c,d,e,f). 

Figure 5: The relative phase-errors. Notation and labelling like in Figure 3. 

Figure 0S The density distribution functions. Notation and labelling like in Figure 4. 
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Figure lb 





Figure lc 












































